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1. Introduction 

Matter-wave vortices represent fundamental nonlinear macroscopic excitations of Bose- 
Einstein condensates (BECs); see e.g. the relevant reviews [Tf2"f3fHl5ll6]rT] • These 
structures are characterized by their nonzero topological charge S, the phase dislocation 
and jump by 2irS induced by the vorticity and the concomitant vanishing of the BEC 
density at the vortex core. Experimental observation of matter-wave vortices was first 
reported in Ref. [8], using a phase-imprinting method between two hyperfine spin states 
of a 87 Rb BEC [9]. Other techniques for the generation of vortices have also been 
studied theoretically and implemented in experiments. In particular, stirring the BEC 
[TO] above a certain critical angular speed [TTfl2|ll3j is an extremely efficient method 
for producing a few vortices [13] or vortex lattices [14]. Other techniques include the 
supercritical dragging of an obstacle through the BEC [T5lll6fl7] . as well as the nonlinear 
interference of condensate fragments [T8|19|20|[2T] . In the above studies, vortices were 
singly- charged i.e., with a topological charge S = 1; higher-charged vortices with 
S > 1 may also be created experimentally [22|23] and could, in principle, be stable 
under appropriate conditions [24f25] . Considerable effort has been dedicated to the 
investigation of the stability of such higher charge structures [26112711281(29] . Nevertheless, 
such higher-charged vortices are typically far less robust than the fundamental 5=1 
vortex that is of interest here. 

In this work we systematically study singly-charged vortices in a two-dimensional 
(2D) — so-called disk-shaped — BEC from a spectral (i.e., Bogoliubov-de Gennes) point 
of view. In particular, first we focus on the so-called anomalous mode of the Bogoliubov 
theory, characterized by negative energy [30] or negative Krein-sign [31] , and elucidate its 
connection with the precessional motion of the vortex, if displaced from its equilibrium 
position i.e., the trap center. Next, we will study how this mode is affected by the 
presence of different kinds of perturbations. The perturbations we consider here arise 
from inhomogeneous interatomic interactions, so-called collisional inhomogeneities, and 
finite-temperature induced dissipation. 

Interatomic interactions, characterized by the s-wave scattering length, become 
spatially (or temporally) varying by employing either magnetic [32f3 3j or optical 
Feshbach resonances [34.35.36j in a very broad range. This remarkable flexibility on the 
manipulation of the effective mean-field nonlinearity of BECs, has inspired a significant 
number of experimental and theoretical studies. Herein, we will focus on the more 
recently proposed "collisionally inhomogeneous" BECs, characterized by a spatially- 
dependent scattering length. In such settings, many interesting phenomena have been 
predicted, including adiabatic compression of matter- waves [371138] . Bloch oscillations 
of solitons [37], soliton emission and atom lasers [39], enhancement of transmittivity 
of matter- waves through barriers [40f4T] . dynamical trapping of solitons [40], stable 
condensates exhibiting both attractive and repulsive interatomic interactions [42] and 
the derealization transition of matter waves [43]. Here we will examine how harmonic 
spatial variations of the scattering length, inducing a sort of a nonlinear optical lattice in 
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the system, affect the stability and ensuing dynamics of the vortex. Interestingly, we find 
that the anomalous mode of the vortex (located at the origin) is differently affected by 
cosinusoidally (the vortex is located at a maximum of the nonlinearity) and sinusoidally 
(the vortex is located at a local minimum of the nonlinearity) varying nonlinearities i.e., 
the phase of the nonlinearity's spatial variation at the vortex location plays a crucial 
role in the ensuing stability properties. This turns out to be the most critical element 
of influence within this setting. In the former case the vortex is stable, while in the 
latter the vortex is subject to an oscillatory instability, emerging by the collision of the 
anomalous mode with another eigenmode of the system. 

We also consider in our study the effect of dissipative perturbations on the vortex 
dynamics motivated by considerations of the coherent structure's interaction with the 
thermal cloud. Here we will adopt a simple phenomenological model relying on the 
inclusion of a phenomenological damping in the mean-field model, first introduced by 
Pitaevskii [11] and subsequently used in various works to describe, e.g., decoherence 
[15] and growth [16] of BECs, damping of collective excitations of BECs [17] . vortex 
lattice growth jl8fl9] . vortex dynamics [50] (see also [T9|20] ). and decay of dark solitons 
[51] . Importantly, inclusion of such a phenomenological damping in the Gross-Pitaevskii 
equation (GPE) can be justified from a microscopic perspective (see, e.g., the recent 
review [52]). Herein, we will show how such a finite-temperature motivated dissipation 
affects the statics and dynamics of the vortex, by leading its anomalous mode to 
become immediately unstable. Despite the relatively simple and phenomenological 
nature of the model, we will see that its results will bear significant similarities to the 
phenomenology of more complex dynamical models of the relevant interactions, allowing 
us to understand qualitatively the origin of the observed dynamical features. We will 
also present some interesting twists that may arise when the combined effect of thermal 
dissipation and spatially-dependent interatomic interactions is considered. 

The paper is structured as follows. In Section [2] we present our analytical 
considerations in connection to the spectrum of a S = 1 vortex and its precession 
frequency in the trap. In Section [3l we examine numerically the validity of the 
analytical predictions, but also how these are modified in the presence of additional 
perturbations such as the spatially dependent nonlinearity, or the finite-temperature 
induced dissipative perturbation. This is done both through a systematic analysis of 
the Bogoliubov-de Gennes equations, as well as through the direct numerical simulations 
of the pertinent GPE models. Finally, in SectionH] we summarize our results and discuss 
directions for future studies. 

2. Analytical Results 

We first consider the simplest case in our study, namely a matter-wave vortex in a 
2D BEC confined in a harmonic trap. It is well-known that in this setting the singly 
charged vortex will have precisely one anomalous mode [30]; this mode, characterized 
by a negative energy, is also known as mode of negative Krein sign (or signature) in the 
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mathematical literature [3T]; see also below for an explicit mathematical definition. In 
Ref. [30] (see also the review 0]), it was argued that the single negative energy mode 
with S = 1 (which is of interest here) is responsible for the precessional motion of 
the vortex in the trap (in addition to being relevant for other processes such as vortex 
nucleation) . 

One of the key purposes in our study is to consider the precession in the setting 
described above. The three-dimensional (3D) analogue of this setting has been 
considered and studied analytically by means of the matched asymptotics technique 
in Ref. [53], while the 2D case has been studied by means of a variational approach in 
Ref. [51] (see more details below). Here, employing the matched asymptotics method, we 
derive an expression for the precession frequency in the 2D case, and provide a detailed 
comparison of this result with numerics pertaining to the study of the anomalous mode. 

The model under consideration is the (2 + l)-dimensional GPE [7], 

h 2 

ihdfU = A'u + V(r')u + g2DW\ 2 u — u'u. (1) 

2m 

Here, u(x',y',t') is the macroscopic wave function of the disk-shaped BEC, A' is 
the 2D Laplacian, r' = a/ x' 2 + y' 2 is the radial variable, V(r') = \u 2 r' 2 is the 
harmonic trapping potential in the in-plane direction, // the chemical potential and 
92D = 93d/^Oj Z = 2y/2TTaa z fvuj z is the effectively 2D nonlinear coefficient where a 
is the scattering length and a z , u z are the transverse (strongly confining) harmonic 
oscillator length and trapping frequency, respectively. Measuring length in units of a z 
and frequencies in units of u z Eq. ([TJ can be expressed in the following dimensionless 
form, 

id t u = — Au + V(r)u + g\u\ 2 u — fiu. (2) 
2 

Here, A is the 2D Laplacian of the rescaled variables, r is the rescaled radial variable, 
V(r) = |f2 2 r 2 is the harmonic trapping potential with fl being measured in units of u z , 
g is the normalized strength of the interatomic interactions (which we set to g — 1 for 
the analytical considerations of this section), and \x is the chemical potential measured 
in units of ftw z . In order to study the effect of the potential on the vortex, we will 
follow the lines of Ref. [55] (see also Ref. [56] for similar work in the context of optics) 
and use a matched asymptotics approach between an inner and an outer perturbative 
solution leading to the following equations of motion (for a more detailed derivation see 
the appendix) for the vortex center (x v ,y v ): 

x v = ^log(^)^ (3) 

*»= -| log KH (4) 

where A is an appropriate numerical factor (detailed comparison with numerics yields 
very good agreement in the Thomas-Fermi regime e.g. for A w 8.88 w 2y/2n, see below). 
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These results suggest a precession of the vortex in the harmonic trap with a frequency 



which, as suggested by the subscript ("an" stands for anomalous), should coincide with 
the eigenfrequency of the anomalous mode of the Bogoliubov spectrum. The anomalous 
mode eigenfrequency can readily be obtained through a standard Bogoliubov-de Gennes 
(BdG) analysis. This analysis involves the derivation of the BdG equations, which stem 
from a linearization of the GPE (CQ) around the vortex solution uq by using the ansatz 



The subsequent solution of the ensuing BdG eigenproblem yields the eigenfunctions 
{a(x,y),b(x,y)} and eigenfrequencies u. 

Note that due to the Hamiltonian nature of the system, if u is an eigenfrequency 
of the Bogoliubov spectrum, so are —u, u* and —u*. Notice that a linearly stable 
configuration is tantamount to Im(oj) = 0, i.e., all eigenfrequencies being real. 

An important quantity resulting from the BdG analysis is the amount of energy 
carried by the normal mode with eigenfrequency u, namely 



The sign of this quantity, known as Krein sign [57], is a topological property of each 
eigenmode. For one of the eigenvalues of each double pair this sign is negative. The 
corresponding mode is called negative energy mode (in the physical literature) [58] or 
mode with negative Krein signature (in the mathematical literature) [57]. Practically, 
this means if it becomes resonant with a mode with positive Krein signature then, 
in most cases, there appear complex frequencies in the excitation spectrum, i.e., a 
dynamical instability occurs [57]. The eigenvalues with negative Krein signature are 
actually associated with the anomalous modes [1] appearing in the BdG spectrum. 

In order to compare our results to the ones obtained in earlier works, we should 
mention that a similar setup was investigated in Ref. [59] for finite displacements of 
the vortex from the center of the trap and in Ref. [M] (by means of a variational 
approximation). In the latter one the frequency of the anomalous mode was derived with 
a similar functional form. However, in Ref. [54], the case of a BEC unbounded in the 
axial direction (oo z = 0) was considered and, as a result, the constant was found to take 
a different value, A = 2. It is also worth noting that the connection between quantum 
fluctuations and anomalous modes of matter-wave vortices under Magnus forces was 
considered in Ref. [60] 

It is important, at this stage, to make a few comments regarding the nature of the 
Bogoliubov spectrum resulting from the linearization around the vortex. The system 
at hand, namely the disk-shaped condensate carrying the vortex, is not in the ground 
state (a similar situation occurs in the ID analogue of the system, namely a quasi-lD 
BEC carrying a dark soliton). The existence of the anomalous mode, characterized 
by negative energy, indicates that the vortex (and the dark soliton in the ID case) is 




(5) 




(6) 




(7) 



Stability and Dynamics of matter-wave vortices 



6 



thermodynamically unstable and, in the presence of dissipation, the system is driven 
towards a lower energy configuration, namely the ground state. Also, the eigenfrequency 
of the anomalous mode of the vortex (similarly to the case of the dark soliton |61j ) 
bifurcates from its value u = Q (the linear oscillation with the trap frequency) in the 
linear limit — where the vortex is represented by the linear superposition |1, 0) + z|0, 1), 
where \m, n) denotes the m-th linear eigenstate of the quantum harmonic oscillator 
along the re-direction and ra-th one along the y-direction (see discussion in Section 
13.11 and Fig. [T]). Generally, the anomalous mode (in both ID and 2D cases) is the 
lowest excitation frequency of the system and the only one below the trap frequency, 
which is associated with the doubly degenerate (in the two-dimensional case) Kohn 
mode corresponding to the dipolar motion of the condensate |62j. However, in the 
case of the vortex (and contrary to what is the case for the dark soliton), there is one 
more frequency, which may be smaller than w an , at least for small chemical potentials. 
However the corresponding mode grows monotonically away from the origin and thus 
becomes larger for increasing chemical potential than the anomalous mode and the Kohn 
mode. This frequency was described through a small parameter expansion in Section 
V.B of Ref. [63]. The relevant eigenfrequency is given (in our units) by the following 
expression, 

u « n - 2Q, (8) 

which becomes increasingly more accurate as fi — > 2Q. This is in contrast to the case 
for the expression of the precession frequency Eq.(jSJ), which should be increasingly more 
accurate in the Thomas-Fermi (TF) limit, corresponding to large /i. 

We now turn to numerical investigations in order to examine the validity of our 
results in the case of the parabolic trap for constant nonlinearity strength, as well as to 
generalize them to settings which are less straightforward to consider by analytical 
means. The results will be partitioned in two subsections: firstly, we will provide 
bifurcation results from the BdG analysis, and subsequently, we will also test the BdG 
predictions against full numerical integration of Eq. (JT]). 

3. Numerical Results 

3.1. Bogoliubov-de Gennes Analysis 

We start with the case of a singly-charged vortex in a harmonically confined BEC with 
homogeneous interatomic interactions (i.e., g = 1). In Fig. [T]we show the numerically 
obtained eigenfrequency u of the Bogoliubov spectrum as (blue) solid lines as a function 
of the chemical potential \x and the analytical predictions of Eq. (J5]) (green) dashed- 
dotted line and Eq. flH]) (red) dashed line. The numerically obtained frequencies are 
real denoting that the system is dynamically stable. We observe that in accordance to 
the analytical predictions, the lowest modes are (i) the one monotonically increasing 
away from zero and (ii) the anomalous mode, connected to the vortex precession (see 
previous section), bifurcating from the constant dipolar mode in the liner limit. For the 
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Figure 1. (Color online) The eigenfrequency cj of the Bogoliubov spectrum for a 
harmonically confined singly-charged vortex as a function of the chemical potential 
/x (for a trap strength f2 = 0.2). Theoretical predictions are given by u> = 
(/x — 20) [(red) dashed line] for the mode monotonically increasing from zero and 
u = (Sl 2 /(2/i)) log(Afi/Sl) [(green) dash-dotted line]; the constant A was chosen to be 
2y/2n « 8.886. 

monotonically increasing mode, we notice that the non-radial nature of the solutions at 
hand (due to their phase profile) leads to the absence of additional symmetries of the 
eigenvalue problem away from the linear limit. The only symmetry generally present is 
that of the phase or gauge invariance, associated with the conservation of the number 
of particles. This sustains a pair of linearization eigenfrequencies at u = 0, but as 
discussed in Ref. [53], at the linear limit the dimension of the corresponding kernel is 
4, hence an eigenfrequency pair should depart from the origin (at least for small /x) 
according to Eq. (jE}. As observed in Fig. [TJ this prediction is in good agreement with 
the numerical results. Naturally, deviations are observed for larger chemical potentials. 
On the other hand, as concerns the precession frequency, we notice its monotonically 
decreasing dependence on /x for given Q (the latter, was set to Q = 0.2 in Fig. [TJ, 
its bifurcation from the Kohn mode eigenfrequency limit and its excellent agreement 
with the theoretical prediction in the TF limit (for all /x > 1). We note in passing 
that in this two-dimensional case, a pair of Kohn modes can be seen to be preserved 
at u = f2 = 0.2, being associated with the dipole oscillations of the condensate along 
the two spatial directions, while the fourth mode at 0.2 in the linear limit results in a 
monotonically growing eigenfrequency, as /x is increased. 

We have also tested the validity of the analytical predictions concerning the two 
lowest eigenfrequencies for different values of Q and as a function of /x (see Fig. [2]). 
Once again, a very good agreement of the two asymptotic theoretical descriptions in 
their respective limits is found. 

We now consider an interesting modification to this picture, arising from a 
consideration of inhomogeneous interatomic interactions, described by a spatially- 
dependent scattering length a(x,y) (see, e.g., the recent special volume [64]). Here, we 
will consider the effect of a periodic variation of the nonlinearity strength, g = g(x,y) 
(i.e., a sort of a nonlinear optical lattice) on the spectrum of a vortex. We will also draw 
parallels with similar spectral implications in the setting of a linear periodic potential 
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Figure 2. (Color online) Top: dependence of the anomalous mode eigenfrequency on 
the chemical potential for different trap frequencies f2. The data points are interpolated 
by the functions ,fn{p) = (fi 2 /(2/i)) log(2v / 27r/i/J7). Bottom: similar to top panel, 
but for the mode bifurcating from zero, as compared to the theoretical prediction 
gn(fJ>) = V - 20. 



analyzed in Ref. [65] . 

In Fig. [31 we study the case of a cosinusoidal variation of the nonlinearity strength, 
namely, g(x,y) = 1 + s (cos 2 (irx/4) + cos 2 (it y / 4)), monitoring the vortex spectrum as 
a function of the chemical potential, where s is the strength of the oscillation. In the 
same figure, the typical form of the density and phase of the wave function (the former 
showcasing spatial variation dictated by the corresponding variation of the scattering 
length, and the latter demonstrating the vortex structure of the configuration), as well 
as the Bogoliubov excitation spectrum, are also illustrated. We notice that while most 
of the relevant eigenfrequencies are only very weakly affected by the spatially-dependent 
nonlinearity, the one which is most dramatically affected is that of the anomalous mode. 
The comparison of the s = case of Fig. [T] (blue solid lines) with the red dashed line 
of s = 0.1 and the green dash-dotted of s = 0.3 illustrates that the latter two not only 
approach zero, but rather cross it at a finite value of /i. For s = 0.1, the anomalous mode 
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hits the origin of the spectral plane at fi — 2.61, while for s = 0.3 at [i = 1.56. However, 
it is perhaps even more remarkable that this collision does not produce an instability 
through an imaginary eigenfrequency (real eigenvalue) pair, but rather maintains the 
stability of the configuration (the eigenfrequencies appear to go through each other). 
Generally, it can be seen that the trend of increasing the oscillation strength in the 
cosinusoidal case leads to a more rapid decrease of the anomalous mode eigenfrequency 
with fi and an "earlier" collision (i.e., occurring for smaller /i) with the spectral plane 
origin. The present study focuses on a periodicity (wavelength) of the nonlinearity that 
is larger than the size (core) of the vortex. All throughout this long wavelength regime 
the spectral results are qualitatively the same. The case pertaining to wavelengths of the 
nonlinearity comparable or smaller than the size of the vortex falls outside of the scope of 
the current manuscript and will be studied further in a future work. Nonetheless, it can 
be anticipated that for small enough wavelengths compared to the core of the vortex, the 
spatial modulation of the nonlinearity will effectively, through spatial homogenization, 
act as a constant nonlinearity (possibly shifted from its original g = 1 value) in a manner 
akin to the effects of (linear) periodic potentials generated by optical lattices acting on 
harmonically trapped dark solitons [66] . 

It is now interesting to turn to the case of the sinusoidal modulation of the 
nonlinearity strength, namely g(x,y) — 1 + s (sin 2 (7rx/4) + sin 2 (7n//4)). In this case, 
as observed in Fig. HI the fundamental difference is that the anomalous mode is larger 
than that of the homogeneous interactions case (g = 1). More importantly perhaps, its 
dependence can also be non-monotonic, resulting in the increase of the corresponding 
eigenfrequency for a chemical potential \i > 1. Consequently, this raises the possibility 
of collision of the relevant eigenmode with other modes bifurcating from u = Q for 
sufficiently large /i (see [green] cross fora fi ~ 2.11 in the case of s = 0.3 considered 
in the figure). This, in turn, produces an instability due to the opposite Krein sign of 
the colliding modes, yielding a quartet of complex eigenfrequencies. The (positive) 
imaginary part of the latter, is shown in the bottom panel of Fig. HI see also the 
second and third row of panels for a typical profile and spectral plane of the relevant 
configuration. 

The above features of the anomalous mode seem structurally similar to the linear 
periodic potential case, where again the cosinusoidal case was found to be dynamically 
stable, while the sinusoidal one was unstable beyond a critical lattice strength [65l6|67j . 




These results also motivate an investigation of how this phenomenology may be modified 
in the presence of a dissipative term. 

In this case, the pertinent model is the so-called dissipative GPE, which can be 
expressed in the following dimensionless form: 



where the dimensionless parameter 7 can be associated with the system's temperature 





(9) 
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Figure 3. (Color online) The case of a cosinusoidal variation of the nonlinearity 
strength for a trap strength il = 0.2. The top panel is similar to Fig. [T] showing the 
cigenfrequency uj of the Bogoliubov spectrum as a function of the chemical potential /i. 
Comparison of g(x,y) = 1 + s (cos 2 (irx / 4) + cos 2 (iry / 4)) , with s = 0.1 [(red) dashed 
line] and s = 0.3 [(green) dash-dotted line], with the case of s = [(blue) solid line). 
The middle panels show contour plots of the density (left) and phase (right) of the wave 
function, while the bottom panel shows the respective Bogoliubov excitation spectrum 
(real vs. imaginary part of the eigenfrequency w, where instability would correspond 
to the existence of eigcnfrequcncics with Im(cj) ^ 0). The chemical potential is fi = 4 
(approaching the Thomas- Fermi limit). 



in SI units according to p£6|H9] (see also [52] ) 

Ama 2 k B T 

7 = Gx « — ' 10 

7Th Z 

with ks being Boltzmann's constant and the heuristically introduced dimensionless 
prefactor G ~ 3. Note that in the dissipative model the interaction between the 
thermal cloud and the condensate is only modeled by particle exchange resulting in 
the dissipative factor 7; we should note that physically, the relevant case is that of 
7 < 1, although for illustration purposes we will occasionally show also the results 
away from that regime. The chemical potential and trap strength in Eq. ((HD are set to 
the values \i = 1 and Q = 0.2 (per the above discussion, it is understood how different 
jj, and Q will modify the relevant phenomenology). 

In Fig. [5l we show the BdG spectrum of a vortex for a case of 7 = (zero 
temperature, i.e., no dissipation) and for the case of 7 = 0.2 (finite temperature, i.e., 
dissipation). It is clear that the lowest frequency mode of the condensate (which for 
fi = 1 is the anomalous mode) is the one that, for nonzero values of 7, immediately 
acquires a positive imaginary eigenfrequency, contrary to what is the case for all other 
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Figure 4. (Color online) Similar to Fig. |3l but for the case of a sinusoidal variation 
of the nonlinearity strength. The first three rows of panels are analogous to the panels 
of Fig. H The case of g(x, y) = 1 + s (sin 2 (7rx/4) + sin 2 (7n//4)) , for s = 0.3 [(red) 
dashed line] is compared to that of s = [(blue) solid line]. The bottom panel shows 
the imaginary part of the complex eigenfrequency; the oscillatory instability arises for 
y > 2.11. 
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eigenmodes of the system. In fact, precisely this property was rigorously proved in 
Ref. [31] for negative Krein sign eigenmodes, namely that their bifurcation upon such 
dissipative perturbations happens oppositely to that of all other modes (of positive 
energy) of the system. This remarkable feature is directly consonant with the property 
of this excited state of the system resulting, via the effect of dissipation (and through 
the complex nature of the relevant eigenmode) eventually into the ground state of the 
system. Moreover, notice that this complex eigenmode also implies the combination 
of growing amplitude with the previously analyzed precessional motion, leading to 
the spiraling of the vortex core toward the edges of the TF cloud and its eventual 
disappearance, in favor of the ground state of the system (see below). 

Lastly, let us investigate the effect of a periodic modulation of the nonlinearity on 
the stability of the system for the finite-temperature case. For a periodic cosinusoidal 
modulation of the nonlinearity the imaginary part of all eigenfrequencies vanishes and 
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Figure 5. (Color online) The top panel shows the immediate acquisition of a non- 
vanishing imaginary part of the eigenfrcquency associated with the anomalous mode, 
as soon as the temperature-dependent dissipative prcfactor 7 becomes nonzero. The 
bottom panel highlights the special behavior of the anomalous mode by illustrating the 
BdG spectrum for the cases of 7 = (red crosses) and that of 7 = 0.2 (blue circles). 
Notice how in the latter case the eigenfrequencies form nearly two symmetric arcs in 
the negative imaginary half-plane. 



the system is stable for 7 = as discussed above. The top panel of Fig. [6] shows the 
maximal imaginary part of the eigenfrequency as a function of 7 for different chemical 
potentials \x for g(x,y) = 1 + s (cos 2 (7rx/4) + cos 2 (-7n//4)), with s = 0.3. For 7 = 
the imaginary part of the eigenfrequency vanishes for all cases. However, for small 
/j, the imaginary part of the eigenfrequency becomes non-zero immediately, similar to 
the case of a constant nonlinearity strength. On the other hand, for large chemical 
potential the maximal imaginary part of the eigenfrequencies remains zero independent 
of 7. Thus, the system remains stable even in the presence of dissipation. This behavior 
can be understood by investigation of Fig. HJ The occurrence of a positive imaginary 
part of the eigenfrequencies is due to the fact that the anomalous mode is of negative 
Krein sign. However, for the case of a cosinusoidal modulation of the nonlinearity 
one observes that the value of the frequency of the anomalous mode decreases with 
increasing /i and, finally, even crosses the origin. At that critical point, the frequency 
curve shows a crossover with its opposite- value companion (of the same pair). However 
the latter mode has a positive Krein sign and therefore (since all negative signatures 
arise for negative frequencies, and positive signatures for positive frequencies), the 
imaginary parts of the eigenfrequencies become negative in the case of nonzero 7. The 
bottom panel in Fig. [6] provides an overview of the eigenfrequencies for 7 = 0.2 and 
/i = 1.6. All eigenfrequencies were shifted further into the negative imaginary half-plane 
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Figure 6. (Color online) The top panel shows the immediate acquisition of a non 
vanishing imaginary part of the eigenfrequency associated with the anomalous mode, as 
soon as the temperature-dependent dissipative prefactor 7 is nonzero for small chemical 
potentials, but no acquisition of an imaginary part for large chemical potentials. The 
bottom panel gives an overview of the eigenfrequencies for 7 = 0.2 and /i = 1.6. 
The maximum imaginary part is zero. The inset shows that the eigenfrequencies 
corresponding to the anomalous mode got shifted into the negative imaginary half- 
plane. 



in comparison to the corresponding panel in Fig. |5j Importantly, the eigenfrequencies 
corresponding to the (formerly) anomalous mode got shifted into the negative imaginary 
half-plane as is shown in the inset. 

3.2. Direct Numerical Simulations 

In this section we show results obtained by direct numerical integration of Eq. (pQ) 
starting with different initial states containing a single vortex. In order to determine 
the position of the vortex as a function of time we first compute the fluid velocity [16] 

i u*Vu — uS7u* ,„ „. 



The fluid vorticity is then defined as oo vor = V x v s . Due to our setup, the direction of 
the fluid vorticity is always the z-direction and, therefore, we can treat this quantity as 
a scalar. Furthermore, we investigate single vortex states leading to a single maximum 
of the fluid vorticity at the position of the vortex. This allows us to determine the 
position of the vortex by determining the center of mass of the vorticity u VOT . 
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Figure 7. (Color online) The top panel shows the trajectory of the vortex for g = 1 
and \i = 3 obtained by direct integration (crosses) and the theoretical prediction (solid 
line) obtained by solving Eqs. (J3J (|4j) . Notice the excellent agreement between the 
two. The bottom panel shows the corresponding trajectory for the case taking into 
account dissipation, namely integrating Eq. ([9]) with 7 = 0.2. 

Figure [7] shows the evolution of a single vortex for g — 1. We displaced the vortex 
initially from the center of the trap to (xo,yo) = (—1.5, 0) and propagated the state 
numerically using Eq. ([[]). The thus obtained results are compared to the solutions of 
Eqs. ©-(U x = x cos(Ct) and y = y sin(Ct) with C = (tt 2 /(2/j,)) log(A///fi) and the 
initial position (xo,yo). The theoretical predictions agree very well with our numerical 
findings: the vortex oscillates around the center of the trap with constant frequency 
and radius (see top panel). The bottom panel shows the trajectory for the case of 
constant nonlinearity g = 1 but for finite temperature (i.e., including dissipation). The 
results shown were obtained by direct numerical integration of Eq. (Q for 7 = 0.2, with 
the initial condition being a slightly perturbed eigenstate of the system. Due to the 
instability of the system this small perturbation leads to the spiraling out of the vortex, 
as is physically anticipated in the presence of finite temperature [68J ; we note in passing 
that this work contains a detailed model from microscopic first principles that illustrates 
a similar phenomenology upon a spatially dependent inclusion of the coupling of the 
condensate with the thermal cloud. 
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Figure 8. (Color online) The trajectory of a vortex for y = 3 for the case 
g(x, y) = 1 + s (sin 2 (7ra:/4) + sin 2 (7rj//4)) with s = 0.3. In the top panel the trajectory 
is plotted on top of the profile of the coupling g(x, y), whereas the bottom panel shows 
the trajectory as a function of time. 

Figure E] shows the trajectory of a vortex for the case of a periodically 
modulated sinusoidal nonlinearity, g(x,y) = 1 + s (sin 2 (7nr/4) + sin 2 (7n//4)). The initial 
configuration is a slightly perturbed eigenstate leading to a small shift of the position of 
the vortex. Due to the instability of the sinusoidal g{x,y) landscape, the vortex spirals 
outwards initially, but then spirals inwards after reaching a region with approximately 
constant nonlinearity. Subsequently the vortex follows a series of such alternating 
(spiraling first outwards, and then inwards) cycles in an apparently quasi-periodic orbit. 

Figure [9] shows the trajectory of a vortex for the case of a periodically modulated 
cosinusoidal nonlinearity, g(x,y) = 1 + s (cos 2 (7rx/4) + cos 2 (iry/ '4)), without dissipation 
(left panel) and with dissipation (right panel). In this case, small perturbations do 
not get amplified since the system is stable. However, a macroscopic displacement of 
the vortex to (xo,yo) = (—1-5,0) leads to the trajectories shown in the figure (see left 
panel). In the case without dissipation the vortex moves outwards (reaching a region 
outside the "square" of the first minima of the nonlinearity) and oscillates around the 
center on a trajectory with roughly constant nonlinearity. For the case with cosinusoidal 
nonlinearity and dissipation (see right panel), the vortex remains stable against small 
perturbations and does not spiral out, contrary to the case of a constant nonlinearity. 
Even for a macroscopical displacement the vortex moves back to the center of the trap 
and remains stable there. This behavior is possible because the effective potential due 
to the spatial variation of the nonlinearity creates the possibility for a metastable vortex 
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Figure 9. (Color online) The trajectory of the vortex for y = 3 for the case 
g{x,y) = 1 + s (cos 2 (-7ra/4) + cos 2 (7ry/4)) with s = 0.3 without dissipation (left) and 
with dissipation (right) with 7 = 0.2. The trajectories are plotted on top of the 
profile of the coupling g(x, y). Notice the fast inward vortex motion in the presence of 
dissipation. 



state even in the presence of dissipation. 

In conclusion, the modulation of the nonlinearity opens up the possibility to 
stabilize the vortex against excitations due to finite temperature effects. This can be 
extremely useful for setups which require stable vortex states for a long period of time 
as, e.g., in the recent work of Ref. [69] which suggests the use of a superposition of two 
counter rotating BECs as a gyroscope. 

4. Conclusions 

In summary, in the present work, we examined the role of anomalous modes in the 
motion of vortices in harmonically confined condensates. We have also focused on the 
settings of spatially dependent scattering lengths and of finite temperature (as well 
as the combination thereof). We found a number of interesting results, including an 
explicit semi-analytical expression for the precession frequency in the trap (by means 
of the matched asymptotics technique), which was found to be in excellent agreement 
with both bifurcation and direct numerical integration results, for different chemical 
potentials and trap frequencies within the Thomas-Fermi regime. 

We subsequently examined how the spectrum (more generally — and the anomalous 
mode in particular) are affected by the presence of spatially-dependent (harmonic) 
interatomic interactions. We found that the latter may induce or avoid instabilities 
depending on the curvature and the strength of the nonlinearity variation. The effect 
of temperature was examined in a simple phenomenological setting which, however, 
still enabled us to observe the thermal instability of the vortex and its rapid spiraling 
towards the edges of the condensate cloud. Intriguingly enough, we also demonstrated 
that the effect of spatially dependent nonlinearities may avoid the thermal instability of 
the vortex by creating a local metastable effective energy minimum wherein the vortex 
can spiral inwards towards the center of the harmonic trap. 

It would be particularly interesting to try to extend both the analytical and the 
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numerical considerations herein towards different directions. On the one hand, it 
is appealing to find similar particle-like equations for the motion (and interaction) 
of multiple vortices within the parabolic trap. On the other hand, it would be 
especially relevant to consider such multi-core realizations in the presence of the thermal 
and spatially dependent nonlinear effects. Yet another direction could be to extend 
considerations presented herein to the case of vortex rings and their dynamics. Such 
studies are presently in progress. 
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Appendix: Equations of Motion for the Vortex 

In this Appendix we detail the effects of the potential on the position of the vortex via 
a matched asymptotics approach between an inner and an outer perturbative solution. 
The inner solution is of the form: 

u(r, 6) = [u (r) + e X (r) cos(0)] ^+^(0^(0)^ (12) 

where e is a formal small parameter associated with the slow speed of precession and 
uo(r) is the radial vortex profile, while x an d V are functions of r, whose asymptotics 
have been detailed in Refs. [551156] (see also Ref. [53J). The outer perturbative solution 
can be obtained by a lowest-order equation for the phase, resulting from a rescaling of 
space, r — > er, and time, t — > e 2 t, namely: 

A6 + F ■ V0 = 0, (13) 

where F = Vlog(|uh| 2 ) (hereafter, boldface is used to denote vectors) and \ub\ 2 is 
the (background, hence the relevant subscript) BEC density in the absence of the 
vortex; notice that the density can be approximated in the Thomas-Fermi (TF) limit 
as \ub\ 2 = n — V(r). 

Interestingly, the similarity of Eq. ( |T3|) to Eq. (20) of Ref. [56] could lead to the 
impression that the detailed formalism of Ref. [56] could be blindly followed giving rise 
to the precessional motion of Eq. (27) therein. However, this turns out to be incorrect. 
Particularly, in Ref. [56] . it was non-generically assumed that F of the outer expansion 
can be accurately approximated by a constant. In our case where F m — f2 2 r//i (for small 
and intermediate distances where the matching with the inner expansion is performed), 
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this approximation is clearly not an appropriate one. Instead, we follow the original 
formulation of Ref. [55], which employs the change of variables <f>(x, y) —> 6(x, y): 

9 X = -S U y - <f>—y\ , (14) 

y = sU x -<f>—x\, (15) 

(we will suppress the iS-dependence hereafter, focusing on singly- charged vortices). 
Then, the equation for reads: 

n 2 n 2 

A<p-—{x<f) x + y<P y )-2—<P = 0, (16) 
which, upon using the transformation = H(r)/ \J [i — V(r), yields 

n 2 

AH H = 27r5(r-r ), (17) 

H 

assuming a point vortex source at r . This leads to the asymptotic behavior H = 
— Ko(m\r — ro | ) , where Kq is the modified Bessel function and m = Vtj \fjl- This should 
be directly compared with Eq. (23) of Ref. [56], showcasing that instead of F r /2 in 
the latter equation, here we have the constant factor m multiplying the distance from 
the vortex core. Once this critical modification is made, the rest of the calculation of 
Ref. [56] can be followed directly, yielding the final result (in the presence of the trap): 

*■= | log KK (18) 
*" = -i log «H (19) 

where the pair (x v ,y v ) defines the location of the vortex center, A is an appropriate 
numerical factor (detailed comparison with numerics yields very good agreement in the 
TF regime e.g. for A Rj 8.88 « 2v^tt, see Sec. II). 

It should be noted that this equation is valid for small displacements from the 
trap center which lends further support to the connection of this dynamics with the 
relevant mode of the BdG analysis. For larger displacements from the vortex center, 
this expression should be appropriately corrected [59] . 
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